# Lerner, Michael. "Local power: Understanding the adoption and design of county wind energy regulation." Review of Policy Research. Forthcoming.
#
# Figure generation file
#
# Contents -------------------------------------
# 1. Map of wind regulation
# 2. Map of years of wind development viability 
# 3. Histogram of setback stringency
# 4. Descriptive statistics table
# 5. Estimated margins plot
# 6. Sensitivity analysis for counties without websites
# 7. Wind turbine construction plot

library(tidyverse)
library(haven)
library(sf)
library(maps)
library(units)
library(pbapply)
library(ggplot2)
library(grid)
library(gtable)
library(spdep)
library(xtable)

setwd("~/Google Drive/My Drive/Research/wind_reg/ROPR Replication Materials/")

load("plot_data_3Aug2021")

# 1. Map of wind regulation -------------------------------------
# Make map for continental US
data('us_states')
states <- st_transform(us_states,st_crs(setbacks))

states$fill <- NA
states$fill[which(tolower(states$NAME) %in% c("colorado", "connecticut","delaware", "kentucky",
                                              "maine", "massachusetts", "minnesota", "new hampshire",
                                              "new mexico",  "new york",  "north dakota",  "ohio",
                                              "oregon",  "rhode island", "south carolina", "south dakota",  "tennessee",
                                              "vermont",  "west virginia",  "wisconsin", "wyoming"))] <- "State siting"
states$fill[which(tolower(states$NAME) %in% c("michigan","new jersey","oklahoma","pennsylvania","texas"))] <- "Municipal siting"

setbacks$setback <- ordered(setbacks$setback,
                            levels = c("none","low","medium","high","ban"),
                            labels = c("None","Low","Medium","High","Ban |"))

setbacks <- st_intersection(setbacks, st_union(states))

plotdat1 <- states %>% filter(is.na(fill)==F) %>% dplyr::select(fill)
plotdat2 <- setbacks %>% mutate(fill = setback) %>% dplyr::select(fill)

plotdat <- rbind(plotdat1,plotdat2)

plotdat$fill <- factor(plotdat$fill,
                       levels = rev(c("State siting","Municipal siting",
                                      "Ban |", "High", "Medium",
                                      "Low", "None")))

colnames(plotdat)[1] <- "Setback\nstringency  "

setbackplot <- ggplot() +
  geom_sf(data = states, color = "black",lwd=0.5,fill=NA) +
  geom_sf(data = plotdat, aes(fill=`Setback\nstringency  `), color = "black",lwd=0.25, alpha = 0.5) +
  scale_fill_manual(values = c(viridis::viridis(n = 5,option = "C"),"#8C510A","#959595"),
                    guide = guide_legend(nrow = 1), na.translate = FALSE) +
  theme(legend.text = element_text(size = 16),
        legend.title = element_text(size = 16, face = "bold"),
        legend.position = "bottom",
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank())

# Make map for continental US
data('alaska')
alaska$fill <- "State siting"

alaska <- alaska %>% dplyr::select(fill)

alaska$fill <- factor(alaska$fill,
                      levels = rev(c("State siting","Municipal siting",
                                     "Ban |", "High", "Medium",
                                     "Low", "None")))

colnames(alaska)[1] <- "Setback\nstringency  "

alaskaplot <- ggplot() +
  geom_sf(data = alaska, color = "black",lwd=0.5,fill=NA) +
  geom_sf(data = alaska, aes(fill=`Setback\nstringency  `), color = "black",lwd=0.25, alpha = 0.5) +
  scale_fill_manual(values = c("#959595")) +
  theme(legend.position = "none",
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank(),
        panel.border = element_rect(color = "black", fill = NA, size = 1))

# Make map for Hawaii
data('hawaii')

load("plot_data_3Aug2021")

hawaii <- st_transform(hawaii,st_crs(setbacks))

hawaii$fill <- NA

setbacks$setback <- ordered(setbacks$setback,
                            levels = c("none","low","medium","high","ban"),
                            labels = c("None","Low","Medium","High","Ban |"))

setbacks <- st_intersection(setbacks, st_union(hawaii))

plotdat1 <- hawaii %>% filter(is.na(fill)==F) %>% dplyr::select(fill)
plotdat2 <- setbacks %>% mutate(fill = setback) %>% dplyr::select(fill)

plotdat <- rbind(plotdat1,plotdat2)

plotdat$fill <- factor(plotdat$fill,
                       levels = rev(c("State siting","Municipal siting",
                                      "Ban |", "High", "Medium",
                                      "Low", "None")))

colnames(plotdat)[1] <- "Setback\nstringency  "

hawaiiplot <- ggplot() +
  geom_sf(data = states %>% filter(NAME == "Hawaii"), color = "black",lwd=0.5,fill=NA) +
  geom_sf(data = plotdat, aes(fill=`Setback\nstringency  `), color = "black",lwd=0.25, alpha = 0.5) +
  scale_fill_manual(values = c(viridis::viridis(n = 5,option = "C"),"#8C510A","#959595"),
                    guide = guide_legend(nrow = 1), na.translate = FALSE) +
  theme(legend.text = element_text(size = 16),
        legend.title = element_text(size = 16, face = "bold"),
        legend.position = "none",
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank(),
        panel.border = element_rect(color = "black", fill = NA, size = 1))

setbackplot <- setbackplot + 
  annotation_custom(ggplotGrob(alaskaplot),
                    xmin = -126, ymin = 5,
                    xmax = -111) +
  annotation_custom(ggplotGrob(hawaiiplot),
                    xmin = -111, ymin = 2,
                    xmax = -104)

ggsave(plot = setbackplot,filename = "Plots/setbackmap.png",width = 12,height = 8,units="in",dpi = 600,bg = "white")

# Calculate % of counties with defined setback of all counties with wind energy ordinance
any_ordinace <- which(is.na(as.numeric(setbacks$setback_detail))==F | setbacks$setback_detail=="ban" | setbacks$setback_detail=="not_specified")
defined_setback <- which(is.na(as.numeric(setbacks$setback_detail))==F | setbacks$setback_detail=="ban")

length(defined_setback)/length(any_ordinace)

# Calculate average stringency of setbacks (1 = none, 5 = ban)
setbacks %>% 
  group_by(state) %>% 
  summarize(m=mean(setback_numeric,na.rm = T)) %>%
  arrange(m)

# 2. Map of years of wind development viability -------------------------------------
load("plot_data_3Aug2021")

## Create years_viable variable for every county
# Load county geometries
county <- st_read("Data/County shapefiles/cb_2018_us_county_20m/",layer = "cb_2018_us_county_20m")
county$fips <- as.numeric(paste0(county$STATEFP,county$COUNTYFP))

# Subset to relevant counties
county <- county %>% 
  filter(STATEFP %in% c("02","15","72")==F) %>% 
  dplyr::select(GEOID, ALAND)

colnames(county)[1:2] <- c("fips","land_area")
county$land_area <- as_units(county$land_area,"m^2")
county$land_area <- set_units(county$land_area,"km^2")
county$fips <- as.numeric(as.character(county$fips))

# Transform to get out of lat-long
county <- county %>% st_transform(3857)

# Load wind potential for 2008, 2014, and "near-future" (2019) technology (80-m turbines)
# Transform CRS and interpolate for each county by summing area of potential land and then dividing by the land area
wp08 <- st_read("Data/Wind potential/2008","pot_wind_cap_080_2008")
wp08 <- st_transform(wp08,crs = st_crs(county))
wp08 <- wp08 %>% select(a30)
wp08 <- st_interpolate_aw(wp08,to = county,extensive = T)
# One county does not interpolate correctly, so fix manually
mismatch_indicator <- NA
counter <- 1
while(all(wp08$geometry[counter]==county$geometry[counter])){
  counter = counter+1
}
wp08$gid <- seq(1,nrow(wp08))
wp08$gid[counter:nrow(wp08)] <- wp08$gid[counter:nrow(wp08)]+1
wp08 <- bind_rows(wp08,data.frame(a30=0, gid = counter))
wp08 <- wp08 %>% arrange(gid)
wp08 <- wp08 %>% select(a30)
colnames(wp08)[1] <- c("wp08_km2")
county$wp08_km2 <- wp08$wp08_km2

wp14 <- st_read("Data/Wind potential/2014","pot_wind_cap_080_current")
wp14 <- st_transform(wp14,crs = st_crs(county))
wp14 <- wp14 %>% select(a30)
wp14 <- st_interpolate_aw(wp14,to = county,extensive = T)
# One county does not interpolate correctly, so fix manually
mismatch_indicator <- NA
counter <- 1
while(all(wp14$geometry[counter]==county$geometry[counter])){
  counter = counter+1
}
wp14$gid <- seq(1,nrow(wp14))
wp14$gid[counter:nrow(wp14)] <- wp14$gid[counter:nrow(wp14)]+1
wp14 <- bind_rows(wp14,data.frame(a30=0, gid = counter))
wp14 <- wp14 %>% arrange(gid)
wp14 <- wp14 %>% select(a30)
colnames(wp14)[1] <- c("wp14_km2")
county$wp14_km2 <- wp14$wp14_km2

wpnf <- st_read("Data/Wind potential/nearfuture","pot_wind_cap_080_near_future")
wpnf <- st_transform(wpnf,crs = st_crs(county))
wpnf <- wpnf %>% select(a30)
wpnf <- st_interpolate_aw(wpnf,to = county,extensive = T)
# One county does not interpolate correctly, so fix manually
mismatch_indicator <- NA
counter <- 1
while(all(wpnf$geometry[counter]==county$geometry[counter])){
  counter = counter+1
}
wpnf$gid <- seq(1,nrow(wpnf))
wpnf$gid[counter:nrow(wpnf)] <- wpnf$gid[counter:nrow(wpnf)]+1
wpnf <- bind_rows(wpnf,data.frame(a30=0, gid = counter))
wpnf <- wpnf %>% arrange(gid)
wpnf <- wpnf %>% select(a30)
colnames(wpnf)[1] <- c("wpnf_km2")
county$wpnf_km2 <- wpnf$wpnf_km2

# Determine duration of time county has had sufficient wind potential for utility-scale project (100)
# How much land does a 100MW, 80m hub height project need? Around 40 km^2. Assuming that 1/4 of landowners say yes, use 160 km^2.
county$wp08_viable <- ifelse(county$wp08_km2 >= 160,1,0)
county$wp14_viable <- ifelse(county$wp14_km2 >= 160,1,0)
county$wpnf_viable <- ifelse(county$wpnf_km2 >= 160,1,0)

county$years_viable <- NA
for(i in 1:nrow(county)){
  if(county$wpnf_viable[i] == 0){
    county$years_viable[i] <- "0 years"
  } else if (county$wpnf_viable[i] == 1 & county$wp14_viable[i] == 0){
    county$years_viable[i] <- "1-5 years"
  } else if (county$wpnf_viable[i] == 1 & county$wp14_viable[i] == 1 & county$wp08_viable[i] == 0){
    county$years_viable[i] <- "6-10 years"
  } else {
    county$years_viable[i] <- ">10 years"
  }
}

county$years_viable <- ordered(county$years_viable,
                               levels = c("0 years", "1-5 years", "6-10 years", ">10 years"))

## Load state shape files
data('us_states')
states <- st_transform(us_states,st_crs(county))

county <- st_intersection(county, st_union(states))

states$siting <- "County siting"
states$siting[which(tolower(states$NAME) %in% c("colorado", "connecticut","delaware", "kentucky",
                                                "maine", "massachusetts", "minnesota", "new hampshire",
                                                "new mexico",  "new york",  "north dakota",  "ohio",
                                                "oregon",  "rhode island",  "south carolina","south dakota",  "tennessee",
                                                "vermont",  "west virginia",  "wisconsin", "wyoming"))] <- "State siting"
states$siting[which(tolower(states$NAME) %in% c("michigan","new jersey","oklahoma","pennsylvania","texas"))] <- "Municipal siting"


colnames(county)[9] <- "Years viable"

viabilityplot <- ggplot() +
  geom_sf(data = county, aes(fill = `Years viable`), color = "black", lwd = 0.25, alpha = 0.5) +
  geom_sf(data = states, color = "grey50",lwd=1,fill=NA) + 
  geom_sf(data = states %>% filter(siting=="County siting"), color = "black",lwd=1,fill=NA) + 
  scale_fill_manual(values = c(viridis::viridis(n = 4,option = "C")),
                    guide = guide_legend(nrow = 1), na.translate = FALSE) +
  theme(legend.text = element_text(size = 16),
        legend.title = element_text(size = 16, face = "bold"),
        legend.position = "bottom",
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank())

ggsave(plot = viabilityplot,filename = "Plots/viabilitymap.png",width = 12,height = 8,units="in",dpi = 600,bg = "white")

# Calculate % currently economically viable land in states with county or municipal siting
local_siting_counties <- which(str_extract(county$fips,pattern = "^[0-9][0-9]") %in% states$GEOID[which(states$siting!="State siting")])
county$local_siting <- 0
county$local_siting[local_siting_counties] <- 1
sum(county$wpnf_km2[which(county$local_siting==1)])/sum(county$wpnf_km2)

# Calculate % currently economically viable land in states with county siting
county_siting_counties <- which(str_extract(county$fips,pattern = "^[0-9][0-9]") %in% states$GEOID[which(states$siting=="County siting")])
county$county_siting <- 0
county$county_siting[county_siting_counties] <- 1
sum(county$wpnf_km2[which(county$county_siting==1)])/sum(county$wpnf_km2)

# Calculate % counties viable 0 years, 1-5 years, 6-10 years, >10 years
table(county$`Years viable`)
table(county$`Years viable`)/nrow(county)

# 3. Histogram of setback stringency -------------------------------------
load("plot_data_3Aug2021")

no_website <- length(which(setbacks$setback_detail=="no_website"))
no_ordinance <- length(which(setbacks$setback_detail=="no_ordinance"))
not_specified <- length(which(setbacks$setback_detail=="not_specified"))
low <- length(which(setbacks$setback=="low"))
medium <- length(which(setbacks$setback=="medium"))
high <- length(which(setbacks$setback=="high"))
ban <- length(which(setbacks$setback_detail=="ban"))
de_facto_ban <- length(which(setbacks$setback=="ban"))-ban

setbackplot <- tibble(setback = c("No website", "No ordinance", "Not specified",
                                  "1ft-399ft","400ft-600ft","601ft-1319ft",">1319ft","Ban"),
                      n_counties = c(no_website, no_ordinance, not_specified, low, medium, high, de_facto_ban, ban))
setbackplot$setback <- factor(setbackplot$setback,levels = c("No website", "No ordinance","Not specified",
                                                             "1ft-399ft","400ft-600ft","601ft-1319ft",">1319ft","Ban")) 

g <- ggplot(setbackplot,aes(x = setback, y = n_counties)) +
  geom_col() + 
  labs(y = "N(counties)", x = NULL) +
  scale_y_continuous(expand = c(0,0), limits = c(0,1200)) +
  geom_label(aes(label = n_counties),nudge_y = 15, size = 3) +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"),
        axis.text = element_text(size = 11))

g <- ggplotGrob(g)

g <- gtable_add_rows(g,unit(30,"mm"))

spaces <- vector(length=16)
spaces[1] <- 0.093
for(i in 2:length(spaces)){
  if(i %in% seq(1,length(spaces),by=2) == F){
    spaces[i] <- spaces[i-1]+0.1
  } else {
    spaces[i] <- spaces[i-1]+0.0105
  }
}


# Save as 8.25x6 PDF
pdf("Plots/setbackhist.pdf", width = 8.25, height = 6, bg = "white", colormodel = "cmyk")

grid.newpage()
grid.draw(g)
grid.draw(textGrob(label = "Setback",
                   hjust = 0.5, y = 0.18, gp = gpar(fontsize = 11)))
grid.draw(textGrob(label = "Stringency",
                   x = 0.1, y = 0.13, gp = gpar(fontsize = 11,fontface="bold")))
grid.polyline(x = c(spaces[c(2,6)],spaces[7:8],spaces[9:10],spaces[11:12],spaces[c(13,16)]),
              y = rep(0.13,10),
              id = sort(rep(1:5,2)),
              gp = gpar(col = viridis::viridis(5,option = "C"), lwd = 5))
grid.draw(textGrob(label = "None\n(n=1169)",
                   x = mean(spaces[c(2,6)]), y = 0.12, vjust = 1, gp = gpar(fontsize = 10)))
grid.draw(textGrob(label = "Low\n(n=60)",
                   x = mean(spaces[c(7,8)]), y = 0.12, vjust = 1, gp = gpar(fontsize = 10)))
grid.draw(textGrob(label = "Medium\n(n=208)",
                   x = mean(spaces[c(9,10)]), y = 0.12, vjust = 1, gp = gpar(fontsize = 10)))
grid.draw(textGrob(label = "High\n(n=28)",
                   x = mean(spaces[c(11:12)]), y = 0.12, vjust = 1,gp = gpar(fontsize = 10)))
grid.draw(textGrob(label = "Ban\n(n=49)",
                   x = mean(spaces[c(13,16)]), y = 0.12, vjust = 1,gp = gpar(fontsize = 10)))

dev.off()

# 4. Descriptive statistics table -------------------------------------
load("plot_data_3Aug2021")

desctable <- as.data.frame(matrix(NA,ncol = 5,nrow = 13))
colnames(desctable) <- c("Variable","Mean","SD","Min","Max")
desctable <- as_tibble(desctable)
desctable$Variable <- c("Setback stringency (ordinal)", 
                        "Years viable", 
                        "Neighbor wind ordinance exists", "Neighbor wind ordinance stringency",
                        "Neighbor wind farm exists", "Neighbor wind farm county (log)",
                        "County GDP (log 2012 USD per capita)", 
                        "County GDP from extraction (log 2012 USD per capita)",
                        "Local Government Jobs (log)",
                        "Population density (log people/km2)", "Republican vote (%, Presidential average 2000-2016)",
                        "Transmission distance (log meters)", "N(counties)")
desctable[,-1] <- apply(desctable[,-1],2,function(x){as.numeric(x)})

desctable[1,-1] <- as.list(c(mean(as.numeric(setbacks$setback),na.rm = T), sd(as.numeric(setbacks$setback),na.rm = T), 
                             min(as.numeric(setbacks$setback),na.rm = T), max(as.numeric(setbacks$setback),na.rm = T)))
desctable[2,-1] <- as.list(c(mean(as.numeric(setbacks$years_viable),na.rm = T),sd(as.numeric(setbacks$years_viable),na.rm = T),
                             min(as.numeric(setbacks$years_viable),na.rm = T),max(as.numeric(setbacks$years_viable),na.rm = T)))
desctable[3,-1] <- as.list(c(mean(as.numeric(setbacks$nb_setback_existence),na.rm = T),sd(as.numeric(setbacks$nb_setback_existence),na.rm = T),
                             min(as.numeric(setbacks$nb_setback_existence),na.rm = T),max(as.numeric(setbacks$nb_setback_existence),na.rm = T)))
desctable[4,-1] <- as.list(c(mean(as.numeric(setbacks$nb_setback_stringency_median),na.rm = T),sd(as.numeric(setbacks$nb_setback_stringency_median),na.rm = T),
                             min(as.numeric(setbacks$nb_setback_stringency_median),na.rm = T),max(as.numeric(setbacks$nb_setback_stringency_median),na.rm = T)))
desctable[5,-1] <- as.list(c(mean(as.numeric(setbacks$nb_wind_farm_existence),na.rm = T),sd(as.numeric(setbacks$nb_wind_farm_existence),na.rm = T),
                             min(as.numeric(setbacks$nb_wind_farm_existence),na.rm = T),max(as.numeric(setbacks$nb_wind_farm_existence),na.rm = T)))
desctable[6,-1] <- as.list(c(mean(as.numeric(setbacks$nb_wind_farm_count),na.rm = T),sd(as.numeric(setbacks$nb_wind_farm_count),na.rm = T),
                             min(as.numeric(setbacks$nb_wind_farm_count),na.rm = T),max(as.numeric(setbacks$nb_wind_farm_count),na.rm = T)))
desctable[7,-1] <- as.list(c(mean(as.numeric(setbacks$county_gdp),na.rm = T),sd(as.numeric(setbacks$county_gdp),na.rm = T),
                             min(as.numeric(setbacks$county_gdp),na.rm = T),max(as.numeric(setbacks$county_gdp),na.rm = T)))
desctable[8,-1] <- as.list(c(mean(as.numeric(setbacks$extraction_gdp),na.rm = T),sd(as.numeric(setbacks$extraction_gdp),na.rm = T),
                             min(as.numeric(setbacks$extraction_gdp),na.rm = T),max(as.numeric(setbacks$extraction_gdp),na.rm = T)))
desctable[9,-1] <- as.list(c(mean(as.numeric(setbacks$local_government_jobs),na.rm = T),sd(as.numeric(setbacks$local_government_jobs),na.rm = T),
                             min(as.numeric(setbacks$local_government_jobs),na.rm = T),max(as.numeric(setbacks$local_government_jobs),na.rm = T)))
desctable[10,-1] <- as.list(c(mean(as.numeric(setbacks$population_density),na.rm = T),sd(as.numeric(setbacks$population_density),na.rm = T),
                              min(as.numeric(setbacks$population_density),na.rm = T),max(as.numeric(setbacks$population_density),na.rm = T)))
desctable[11,-1] <- as.list(c(mean(as.numeric(setbacks$republican_vote),na.rm = T),sd(as.numeric(setbacks$republican_vote),na.rm = T),
                              min(as.numeric(setbacks$republican_vote),na.rm = T),max(as.numeric(setbacks$republican_vote),na.rm = T)))
desctable[12,-1] <- as.list(c(mean(as.numeric(setbacks$transmission_distance),na.rm = T),sd(as.numeric(setbacks$transmission_distance),na.rm = T),
                              min(as.numeric(setbacks$transmission_distance),na.rm = T),max(as.numeric(setbacks$transmission_distance),na.rm = T)))
desctable[13,1] <- paste0("N counties: ",length(unique(setbacks$fips[which(is.na(setbacks$setback)==F)])))

desctable[,-1] <- apply(desctable[,-1],2,function(x){round(x,digits = 3)})

print(xtable(desctable),only.contents = T,file = "Plots/descriptives_table.tex",hline.after = F,timestamp = NULL,comment = F)

# 5. Estimated margins plot -------------------------------------
# Load marginal effects estimates and format
View(read_delim(file = "Results/main_3Aug2021.csv",delim = "&"))


main_margins <- read_delim(file = "Results/main_margins_3Aug2021.csv",delim = "&")
main_margins <- main_margins[c(4:9, 12:13,14:15,18:25,26:27),]
colnames(main_margins) <- str_split(colnames(main_margins),pattern = "\\s{2,}",simplify = T)[,1]
main_margins$margin_oprobit_5 <- str_split(main_margins$margin_oprobit_5,pattern = "\\s{1,}",simplify = T)[,2]
main_margins$margin_oprobit_5[which(main_margins$margin_oprobit_5=="\\\\")] <- NA

# Existence plot
margin_inflate <- main_margins[which(is.na(main_margins$margin_inflate)==F),1:2]
margin_inflate_var <- margin_inflate$X1[seq(1,nrow(margin_inflate),by=2)]
margin_inflate_est <- margin_inflate$margin_inflate[seq(1,nrow(margin_inflate),by=2)]
margin_inflate_se <- margin_inflate$margin_inflate[seq(2,nrow(margin_inflate),by=2)]
margin_inflate_se <- str_remove_all(margin_inflate_se,pattern = "[(]|[)]|‘")

margindat_inflate <- tibble(var = margin_inflate_var, est = margin_inflate_est, se = margin_inflate_se) 

margindat_inflate$est <- as.numeric(str_remove_all(margindat_inflate$est,pattern = "[*]"))
margindat_inflate$se <- as.numeric(margindat_inflate$se)

margindat_inflate$ll <- margindat_inflate$est - (qnorm(0.975)*margindat_inflate$se)
margindat_inflate$ul <- margindat_inflate$est + (qnorm(0.975)*margindat_inflate$se)

margindat_inflate$ll_90 <- margindat_inflate$est - (qnorm(0.95)*margindat_inflate$se)
margindat_inflate$ul_90 <- margindat_inflate$est + (qnorm(0.95)*margindat_inflate$se)

margindat_inflate$var <- c("Years viable\n1-5 years","Years viable\n6-10 years", "Years viable\n>10 years",
                           "Neighbor setback exists", "Neighbor wind farm exists")
margindat_inflate$var <- factor(margindat_inflate$var, 
                                levels = rev(c("Years viable\n1-5 years","Years viable\n6-10 years", "Years viable\n>10 years",
                                               "Neighbor setback exists", "Neighbor wind farm exists")))

existenceplot <- ggplot(margindat_inflate, aes (x = var, y = est*100, ymin = ll*100, ymax = ul*100)) +
  geom_point(size = 2) +
  geom_errorbar(width = 0, lwd = 0.75) +
  geom_errorbar(aes(x = var, y = est*100, ymin = ll_90*100, ymax = ul_90*100),width = 0, lwd = 1.25, color = "black") +
  geom_hline(yintercept = 0, color = "red") +
  scale_y_continuous(limits = c(-10,50),breaks = seq(-10,50,by=10)) +
  labs(y = "% Change in Pr(Policy)",
       x = NULL,
       subtitle = "Policy adoption") +
  coord_flip() + 
  theme_classic() + 
  theme(plot.subtitle = element_text(hjust = 0.5, face = "bold"),
        axis.text = element_text(size = 10))


# Stringency plot
main_margins <- gather(main_margins[,-2],level,value,-X1)
margin_oprobit <- main_margins[-is.na(main_margins$value)==F,]
margin_oprobit_var <- margin_oprobit$X1[seq(1,nrow(margin_oprobit),by=2)]
margin_oprobit_est <- margin_oprobit$value[seq(1,nrow(margin_oprobit),by=2)]
margin_oprobit_se <- margin_oprobit$value[seq(2,nrow(margin_oprobit),by=2)]
margin_oprobit_se <- str_remove_all(margin_oprobit_se,pattern = "[(]|[)]|‘")
margin_oprobit_level <- margin_oprobit$level[seq(1,nrow(margin_oprobit),by=2)]
margin_oprobit_level <- factor(margin_oprobit_level,
                               levels = c("margin_oprobit_5","margin_oprobit_4","margin_oprobit_3","margin_oprobit_2","margin_oprobit_1"),
                               labels = c("Ban","High","Medium","Low","None"))

margindat_oprobit <- tibble(var = margin_oprobit_var, est = margin_oprobit_est, se = margin_oprobit_se, Level = margin_oprobit_level) 

margindat_oprobit$est <- as.numeric(str_remove_all(margindat_oprobit$est,pattern = "[*]"))
margindat_oprobit$se <- as.numeric(margindat_oprobit$se)

margindat_oprobit$ll <- margindat_oprobit$est - (qnorm(0.975)*margindat_oprobit$se)
margindat_oprobit$ul <- margindat_oprobit$est + (qnorm(0.975)*margindat_oprobit$se)

margindat_oprobit$ll_90 <- margindat_oprobit$est - (qnorm(0.95)*margindat_oprobit$se)
margindat_oprobit$ul_90 <- margindat_oprobit$est + (qnorm(0.95)*margindat_oprobit$se)

margindat_oprobit$var <- rep(c("Years viable\n1-5 years","Years viable\n6-10 years","Years viable\n>10 years",
                               "Neighbor setback stringency\nLow","Neighbor setback stringency\nMedium",
                               "Neighbor setback stringency\nHigh","Neighbor setback stringency\nBan",
                               "Neighbor wind farm count"),5)

margindat_oprobit$var <- factor(margindat_oprobit$var,
                                levels = rev(c("Years viable\n1-5 years","Years viable\n6-10 years", "Years viable\n>10 years",
                                               "Neighbor setback stringency\nLow","Neighbor setback stringency\nMedium",
                                               "Neighbor setback stringency\nHigh","Neighbor setback stringency\nBan",
                                               "Neighbor wind farm count")))

stringencyplot <- ggplot(margindat_oprobit, aes (x = var, y = est*100, ymin = ll*100, ymax = ul*100, group = Level, color = Level)) +
  geom_point(position = position_dodge(width = 0.5), size = 2) +
  geom_errorbar(position = position_dodge(width = 0.5),width = 0, lwd = 0.75) +
  geom_errorbar(position = position_dodge(width = 0.5),aes(x = var, y = est*100, ymin = ll_90*100, ymax = ul_90*100),width = 0, lwd = 1.25) +
  geom_hline(yintercept = 0, color = "red") +
  scale_y_continuous(breaks = seq(-40,60,by=10)) +
  scale_color_viridis_d(end = 0.9,direction = -1,name = "Level") +
  guides(color = guide_legend(reverse = T)) +
  labs(y = "% Change in Pr(Level|Policy=1)",
       x = NULL,
       subtitle = "Policy stringency") +
  coord_flip() + 
  theme_classic() +
  theme(plot.subtitle = element_text(hjust = 0.5, face = "bold"),
        legend.position = "right",
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 11),
        axis.text = element_text(size = 10))

g1 <- ggplotGrob(existenceplot)
g2 <- ggplotGrob(stringencyplot)

g <- cbind(g1,g2,size="first")

# Save as 10x8 PDF
pdf("Plots/marginsplot_3Aug21.pdf", width = 10, height = 8, bg = "white", colormodel = "cmyk")

grid.newpage()
grid.draw(g)

dev.off()

# 6. Sensitivity analysis for counties without websites
sensitivitydat <- tibble()
for(x in seq(0,1,by= 0.1)){
  if(x != 0 & x != 1){
    x <- str_remove(as.character(x),pattern = "0")
    temp <- read_dta(paste0("Results/Sensitivity_Fig_D1/bootstrap_250_",x,".dta"))
  } else {
    temp <- read_dta(paste0("Results/Sensitivity_Fig_D1/bootstrap_250_",x,".dta"))
  }
  temp$percent_ban <- as.character(x)
  
  sensitivitydat <- bind_rows(sensitivitydat,temp)
}

sensitivitydat$percent_ban <- as.numeric(sensitivitydat$percent_ban)

sensitivitydat <- sensitivitydat %>%
  na.omit() %>%
  pivot_longer(cols = 2:14,names_to = "var",values_to = "estimate") %>%
  group_by(percent_ban, var) %>%
  summarize(est = median(estimate),
            ll = quantile(estimate,probs = c(0.025)),
            ul = quantile(estimate, probs = c(0.975)),
            ll_90 = quantile(estimate,probs = c(0.05)),
            ul_90 = quantile(estimate, probs = c(0.95)),.groups = 'drop')

sensitivitydat_plot <- ggplot(sensitivitydat %>% filter(var=="b_i_nb_existence") %>% mutate(var = "Neighbor setback exists"), aes(x = percent_ban, y = est)) +
  geom_hline(yintercept = 0, color = "red") +
  facet_wrap(~var) +
  geom_point(position = position_dodge(width = 0.2)) +
  geom_errorbar(aes(ymin = ll, ymax = ul),width = 0,position = position_dodge(width = 0.2)) +
  geom_errorbar(aes(ymin = ll_90, ymax = ul_90),width = 0, lwd = 1, position = position_dodge(width = 0.2)) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1),breaks = seq(0,1,by=0.1)) +
  labs(x = "Percent missing setbacks as wind development bans", y = "Estimate") +
  theme_classic() +
  theme(axis.text = element_text(size = 11),strip.text = element_text(size = 14),axis.title = element_text(size = 14),strip.background = element_blank())

ggsave(sensitivitydat_plot,filename = "Plots/sensitivity_plot.png",width = 8,height = 6,units = "in")


# 7. Wind turbine construction plot -------------------------------------
# First run section 2 of this script to generate years of wind development viability measures
viable_df <- as.data.frame(county) %>% select(fips, `Years viable`) %>% rename(years_viable = `Years viable`)

# Expand for 1982:2019
viable_df <- left_join(
  expand.grid(fips = viable_df$fips,
              year = 1982:2019),
  viable_df)

# Load US Wind Turbine Database
fedw <- st_read("Data/Wind farms/uswtdbSHP/",layer = "uswtdb_v3_0_1_20200514")

# Calculate number of commercial projects (50MW or greater) per county per year
turbine_df <- as.data.frame(fedw) %>%
  filter(p_cap>=50) %>%
  group_by(t_fips, p_year) %>%
  summarize(farm_count = length(unique(p_name)),.groups = 'drop') %>%
  rename(fips = t_fips, year = p_year) %>%
  mutate(fips = as.numeric(fips))

# Join with years viable and calculate number of projects per category per year
turbines_time <- left_join(viable_df,turbine_df) %>%
  group_by(years_viable, year) %>%
  summarize(farm_count = sum(farm_count,na.rm = T)/n()) %>%
  na.omit() %>%
  mutate(years_viable = factor(years_viable,
                               levels = c("0 years","1-5 years","6-10 years",">10 years"),
                               labels = c("0 years viable","1-5 years viable","6-10 years viable",">10 years viable")))

# Make plot
turbine_time_plot <- ggplot(turbines_time , aes(x = year, y = farm_count)) +
  facet_wrap(~years_viable,nrow = 1) +
  geom_col() +
  scale_y_continuous(expand = c(0.0015,0.0015)) +
  labs(x = NULL, y = "Number of projects per county") +
  theme(panel.grid = element_blank(),
        panel.background = element_rect(fill = NA, color = "black", size = 1),panel.spacing = unit(0.75,units = "cm"),
        strip.background = element_rect(fill = NA,colour = "black",size = 1)
  )

ggsave(turbine_time_plot,filename = "Plots/turbine_time_viable.png",width = 8,height = 4,units = "in")


